pkgs <- c("tidyverse", "tidylog", "sp", "raster", "devtools", "RCurl", "sdmTMB", "viridis", "terra", "ncdf4", "chron", "ncdf4", "tidync", "tidyterra")
if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Source code for map plots
devtools::source_url("https://raw.githubusercontent.com/VThunell/Lammska_cod-fr/main/R/functions/map-plot.R")
# Packages not on CRAN
# devtools::install_github("seananderson/ggsidekick") # not on CRAN
library(ggsidekick)
theme_set(theme_sleek())
# Set path
home <- here::here()Make prediction grid for relative prey weight models
Introduction
Make an evenly spaced UTM prediction grid with all spatially varying covariates for the stomach content data and biomass estimate data for the for diet and the biomass data
Read stomach data
d <- read_csv(paste0(home, "/data/stomach/stomachs_v1.csv"))
glimpse(d)Rows: 128,771
Columns: 31
$ pred_ID <chr> "LV_1968_9_19_1_37G4_52_NA_25", "LV_1973_5_11_15_37…
$ herring <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ saduria <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ sprat <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ other <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ other_inverts <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ other_chords <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_tot <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_tot_fb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_sad <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_spr <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_her <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_other <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_other_inverts <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_other_chords <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ year <dbl> 1968, 1973, 1973, 1973, 1973, 1983, 1983, 1971, 197…
$ month <dbl> 9, 5, 5, 5, 5, 12, 12, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6…
$ day <dbl> 19, 11, 11, 11, 11, 16, 16, 19, 19, 19, 19, 19, 19,…
$ day_of_year <dbl> 263, 131, 131, 131, 131, 350, 350, 170, 170, 170, 1…
$ pred_weight <dbl> 130.12439, 38.26520, 38.26520, 45.20873, 32.06226, …
$ pred_weight_fb <dbl> 152.97002, 46.30247, 46.30247, 55.27448, 38.37241, …
$ pred_length <dbl> 25, 17, 17, 18, 16, 60, 77, 31, 33, 33, 30, 32, 28,…
$ pred_weight_source <chr> "estimated", "estimated", "estimated", "estimated",…
$ lat <dbl> 54.25, 54.25, 54.25, 54.25, 54.25, 54.25, 54.25, 54…
$ lon <dbl> 14.5, 15.5, 15.5, 15.5, 15.5, 15.5, 15.5, 14.5, 14.…
$ ices_rect <chr> "37G4", "37G5", "37G5", "37G5", "37G5", "37G5", "37…
$ X <dbl> 467.4220, 532.5780, 532.5780, 532.5780, 532.5780, 5…
$ Y <dbl> 6011.453, 6011.453, 6011.453, 6011.453, 6011.453, 6…
$ data_source <chr> "old_db", "old_db", "old_db", "old_db", "old_db", "…
$ Country <chr> "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV…
$ Depth <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…Make the grid with depth
First make a grid for the biomass data, then subset that based on the extend of the stomach data
x <- d$X
y <- d$Y
z <- chull(x, y)
coords <- cbind(x[z], y[z])
coords <- rbind(coords, coords[1, ]) # close the loop
plot(coords[, 1] ~ coords[, 2]) # plot data
sp_poly <- sp::SpatialPolygons(
list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
)
sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
data = data.frame(ID = 1)
)
cell_width <- 4 # increased from 3 to reduce prediction/index processing time and the NEMO hindcast has a spatial resolution of 3.7
d |>
summarise( n = n(), .by = month) # march is the most common month in the data. Use march for predictions# A tibble: 12 × 2
month n
<dbl> <int>
1 9 7018
2 5 4664
3 12 11496
4 6 9689
5 4 25540
6 3 33919
7 1 10967
8 11 14020
9 10 2627
10 7 2193
11 2 6419
12 8 219pred_grid_tmp <- expand.grid(
X = seq(min(d$X), max(d$X), cell_width),
Y = seq(min(d$Y), max(d$Y), cell_width),
year = c(1963:2022)) |>
mutate(month = 3)
ggplot(pred_grid_tmp |> filter(year == 2019), aes(X, Y)) +
geom_point(size = 0.1) +
theme_void() +
coord_sf()
sp::coordinates(pred_grid_tmp) <- c("X", "Y")
inside <- !is.na(sp::over(pred_grid_tmp, as(sp_poly_df, "SpatialPolygons")))
pred_grid_tmp <- pred_grid_tmp[inside, ]
pred_grid_tmp <- as.data.frame(pred_grid_tmp)
ggplot(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000)) +
geom_point(size = 0.001, alpha = 0.5) +
NULL
plot_map +
geom_point(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
NULLWarning: Removed 5091 rows containing missing values or values outside the scale range
(`geom_point()`).
# Add lat and lon
# Need to go from UTM to lat long for this one...
# https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-r
xy <- as.matrix(pred_grid_tmp |> dplyr::select(X, Y) |> mutate(X = X*1000, Y = Y*1000))
v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84 +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]
pred_grid_tmp$lon <- lonlat[, 1]
pred_grid_tmp$lat <- lonlat[, 2]
ggplot(filter(pred_grid_tmp, year == 1999), aes(lon, lat)) + geom_point()
# Add depth now to remove islands and remaining land
# https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package
# https://emodnet.ec.europa.eu/geoviewer/
dep_raster <- terra::rast(paste0(home, "/data/environment/Mean depth natural colour (with land).nc"))
class(dep_raster)[1] "SpatRaster"
attr(,"package")
[1] "terra"crs(dep_raster, proj = TRUE)[1] "+proj=longlat +datum=WGS84 +no_defs"plot(dep_raster)
pred_grid_tmp$depth <- terra::extract(dep_raster, pred_grid_tmp |> dplyr::select(lon, lat))$elevation
ggplot(pred_grid_tmp, aes(lon, lat, color = depth*-1)) +
geom_point()
pred_grid_tmp$depth <- pred_grid_tmp$depth*-1 # To make depth from negative elevation in relation to sea surface??
pred_grid_tmp <- pred_grid_tmp |> drop_na(depth)
pred_grid_tmp |>
filter(year == 1999) |>
drop_na(depth) |>
#mutate(water = ifelse(depth < 0.00000001, "N", "Y")) |>
ggplot(aes(X*1000, Y*1000, fill = depth)) +
geom_raster() +
NULL
plot_map +
geom_point(data = pred_grid_tmp, aes(X*1000, Y*1000), size = 0.001) +
geom_sf() #Simple feature (https://r-spatial.github.io/sf/articles/sf1.html)Warning: Removed 234060 rows containing missing values or values outside the scale range
(`geom_point()`).
plot_map +
geom_raster(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000, fill = depth)) +
geom_sf()Warning: Removed 3910 rows containing missing values or values outside the scale range
(`geom_raster()`).
Combine oxy, temp, and sal from SMHI hindcast (and models of Copernicus and hindcast) for the pred_grid
# read Hindcast data
hindenv_df <- readRDS(file = paste0(home, "/data/environment/hindcast_1961_2017.rds")) |>
filter(year > 1962) |>
mutate(yearmonth = (year-1963)*12+month)
# We need to model the hindcast for year 2018-2023 and we use march as month.
# make a pred grid for those late observations
predhind <- hindenv_df |>
distinct(lat, lon) |>
replicate_df(time_name = "yearmonth", time_values = seq((2018-1963)*12+3,(2023-1963)*12+3, by=12)) |>
mutate( model = as_factor("hindcast"),
year = floor(1963+(yearmonth/12)),
month = yearmonth %% 12) |> # mod, i.e. the remainder of an integer divide (here month)
add_utm_columns(ll_names = c("lon", "lat"), utm_crs = 32633)
# Predict using the models for oxy, sal and temp. For reference see 1c_combine_oxysaltemp_data.qmd
# load model for oxygen
Mod_oxy <- readRDS(paste0(home, "/R/prepare-data/Mod_oxy.rds"))
# load model for salinity
Mod_sal <- readRDS(paste0(home, "/R/prepare-data/Mod_sal.rds"))
# load model for temperature
Mod_temp <- readRDS(paste0(home, "/R/prepare-data/Mod_temp.rds"))
oxypreds <- predict(Mod_oxy, newdata = predhind) |> mutate(oxy = est) |> dplyr::select(lat, lon, year, month, yearmonth, oxy)
salpreds <- predict(Mod_sal, newdata = predhind) |> mutate(sal = est) |> dplyr::select(lat, lon, year, month, yearmonth, sal)
temppreds <- predict(Mod_temp, newdata = predhind) |> mutate(temp = est) |> dplyr::select(lat, lon, year, month, yearmonth, temp)
# Combine hindcast preds for 2018-2023 with hindcast (1963-2017)
hindcast_allyears <- left_join(salpreds, oxypreds) |> left_join(temppreds) |>
bind_rows(hindenv_df |> dplyr::select(temp, sal, oxy, lat, lon, year, month, yearmonth))Joining with `by = join_by(lat, lon, year, month, yearmonth)`
Joining with `by = join_by(lat, lon, year, month, yearmonth)`map(hindcast_allyears, ~sum(is.na(.))) # No NAs$lat
[1] 0
$lon
[1] 0
$year
[1] 0
$month
[1] 0
$yearmonth
[1] 0
$sal
[1] 0
$oxy
[1] 0
$temp
[1] 0hindcast_allyears |>
summarise(nobs = n(), .by = c(month, year)) |>
arrange(year)# A tibble: 666 × 3
month year nobs
<dbl> <dbl> <int>
1 1 1963 18300
2 2 1963 18300
3 3 1963 18300
4 4 1963 18300
5 5 1963 18300
6 6 1963 18300
7 7 1963 18300
8 8 1963 18300
9 9 1963 18300
10 10 1963 18300
# ℹ 656 more rowshindcast_allyears |>
summarise(mean_oxy = mean(oxy), .by = c(year, month)) |>
ggplot() +
geom_line(aes(year, mean_oxy, color = factor(month)), linetype = "dashed")
hindcast_allyears |>
summarise(mean_temp = mean(temp), .by = c(year, month)) |>
ggplot() +
geom_line(aes(year, mean_temp, color = factor(month)), linetype = "dashed")
Extract oxy, sal and temp from hindcast_allyears to pred_grid
# add yearmonth
pred_grid_tmp2 <- pred_grid_tmp |>
mutate(yearmonth = (year-1963)*12+month)
# function for extraction based on yearmonth
ext_envdat <- function(ayearmonth) {
ext_y = pred_grid_tmp2 |> filter(yearmonth == ayearmonth) |> dplyr::select(lon, lat) # coords from the stomach data, can only be lon and lat for extract()
hindcast_allyears |>
filter(yearmonth == ayearmonth) |>
as_spatraster(xycols = c(2,1), crs = "WGS84", digits = 2) |>
terra::extract(ext_y) |>
dplyr::select(oxy,sal,temp) |>
bind_cols(pred_grid_tmp2 |> filter(yearmonth == ayearmonth))
}
# Use the ext_envdat function to extract env covs values to the observations for hindcast
pred_grid_tmp2 <- unique(pred_grid_tmp2 |> pull(yearmonth)) |>
map(\(x) ext_envdat(x)) |>
list_rbind()
# Many observations lie just outside the hincast_allyears raster extent causing NA values of the covariates
theNAs <- pred_grid_tmp2 |> filter(is.na(oxy) | is.na(sal) | is.na(temp))
theNAs |> dplyr::select(lat,lon) |> distinct() |> summarise(n = n())# 427 latlon combos that are outside the hindcast extent n
1 427# ... and they're along the coast which seems reasonable
plot_map_fc +
geom_point(data = unique(theNAs |> dplyr::select(Y,X)), aes(X*1000, Y*1000)) +
theme_sleek(base_size = 6) +
geom_sf() Warning: Removed 193 rows containing missing values or values outside the scale range
(`geom_point()`).
# the 18300 unique hindcast latlon from which we can use get new lat lon to extract to predgrid with
hc_ull = hindcast_allyears |> filter(yearmonth == 1) |> dplyr::select(lon,lat)
# get the nearest nieghbor
pl <- unique(theNAs |> dplyr::select(lon,lat) ) |>
distance( y = hc_ull, lonlat = TRUE)
# select which row in hc_ull has the smallest distance to each of theNAs
pl2 <- hc_ull[apply(pl, 1, FUN = which.min),] # cf a lengthy and slow dplyr version
# combine the real lat lons with their hc_ull corresponding neighbors, temporarliy rename real lat lon for a left_join below.
pl3 <- bind_cols(unique(theNAs |> dplyr::select(lon,lat) |> rename(reallat = lat,
reallon = lon)), pl2 )
# looks alright
plot_map_fc +
geom_point(data = pl3 |> add_utm_columns(ll_names = c("reallon", "reallat"), utm_crs = 32633), aes(X*1000, Y*1000), color = "blue") +
geom_point(data = pl3 |> add_utm_columns(ll_names = c("lon", "lat"), utm_crs = 32633), aes(X*1000, Y*1000), color = "red", alpha = 0.5) +
theme_sleek(base_size = 6) +
geom_sf() Warning: Removed 193 rows containing missing values or values outside the scale range
(`geom_point()`).Warning: Removed 194 rows containing missing values or values outside the scale range
(`geom_point()`).
# left join in the temporary pl3 latlons
pl4 <- left_join(theNAs |> rename(reallat = lat, reallon = lon), pl3)Joining with `by = join_by(reallon, reallat)`# left join in the covariates based on nearest neighbor lat lon and yearmonth
pl5 <- pl4 |> dplyr::select(-oxy,-sal,-temp,-year) |>
left_join(hindcast_allyears)Joining with `by = join_by(month, yearmonth, lon, lat)`pred_grid_tmp3 <- pl5 |>
dplyr::select(-lat,-lon) |> rename(lat = reallat, lon = reallon) |>
bind_rows(pred_grid_tmp2 |> filter(!is.na(oxy) | !is.na(sal) | !is.na(temp)))
# all are there?(!)
pred_grid_tmp2 |> summarise(n()) == pred_grid_tmp2 |> summarise(n()) n()
[1,] TRUE# no NAs left
pred_grid_tmp3 |> filter(is.na(oxy) | is.na(sal) | is.na(temp)) [1] X Y month lon lat depth yearmonth
[8] year sal oxy temp
<0 rows> (or 0-length row.names)# check plots
plot_map_fc +
geom_tile(data = pred_grid_tmp3, aes(X*1000, Y*1000, fill = oxy), alpha = 0.5) +
theme_sleek(base_size = 6) +
geom_sf() Warning: Removed 234600 rows containing missing values or values outside the scale range
(`geom_tile()`).
plot_map_fc +
geom_tile(data = pred_grid_tmp3, aes(X*1000, Y*1000, fill = temp), alpha = 0.5) +
theme_sleek(base_size = 6) +
geom_sf() Warning: Removed 234600 rows containing missing values or values outside the scale range
(`geom_tile()`).
plot_map_fc +
geom_tile(data = pred_grid_tmp3, aes(X*1000, Y*1000, fill = sal), alpha = 0.5) +
theme_sleek(base_size = 6) +
geom_sf() Warning: Removed 234600 rows containing missing values or values outside the scale range
(`geom_tile()`).
Add cod densities, only for post 1992!
# for the Mod, the fr data is instead used to scale the prediction grid
data_stats <- read_csv(paste0(home, "/data/survey/data_stats.csv")) Rows: 12328 Columns: 43
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): country, haul_id, ices_rect, month_year
dbl (39): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, X,...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.m3 <- readRDS(paste0(home, file = "/data/survey/m3.rds"))
# d_stats is used to scale the prediction grid when predictiong cod densities
data_stats_st <- data_stats |>
summarise(depth_mean = mean(depth, na.omit = TRUE),
depth_sd = sd(depth),
sal_mean = mean(sal, na.omit = TRUE),
sal_sd = sd(sal),
oxy_mean = mean(oxy, na.omit = TRUE),
oxy_sd = sd(oxy),
temp_mean = mean(temp, na.omit = TRUE),
temp_sd = sd(temp))
df_m3 <- pred_grid_tmp3 |>
filter(between(year, 1993, 2021)) |>
drop_na(oxy,
sal,
temp) |>
mutate(depth = ifelse(depth < 0, 0, depth),
depth_sc = (depth - data_stats_st$depth_mean)/data_stats_st$depth_sd,
oxy_sc = (oxy - data_stats_st$oxy_mean)/data_stats_st$oxy_sd,
sal_sc = (sal - data_stats_st$sal_mean)/data_stats_st$sal_sd,
temp_sc = (temp - data_stats_st$temp_mean)/data_stats_st$temp_sd,
depth_sq = depth_sc^2,
temp_sq = temp_sc^2,
year_f = as.factor(year),
quarter_f = as.factor(1))
# Predict cod cpue density with model model from Max
cpue_cod <- predict(m3, newdata = df_m3) |>
mutate(log_density_cod = est) |>
dplyr::select(lat, lon, month, year, log_density_cod)
pred_grid_tmp4 <- pred_grid_tmp3 |>
left_join(cpue_cod, by = c("month", "year", "lon", "lat"))Add ICES areas
# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile(paste0(home, "/data/survey/ICES_StatRec_mapto_ICES_Areas/StatRec_map_Areas_Full_20170124.shp"))
head(shape) OBJECTID ID ICESNAME SOUTH WEST NORTH EAST AREA_KM2 stat_x stat_y Area_27
1 1 47 47A0 59.0 -44 59.5 -43 3178 -43.5 59.25 14.b.2
2 2 48 48A0 59.5 -44 60.0 -43 3132 -43.5 59.75 14.b.2
3 3 49 49A0 60.0 -44 60.5 -43 3085 -43.5 60.25 14.b.2
4 4 50 50A0 60.5 -44 61.0 -43 3038 -43.5 60.75 14.b.2
5 5 51 51A0 61.0 -44 61.5 -43 2991 -43.5 61.25 14.b.2
6 6 52 52A0 61.5 -44 62.0 -43 2944 -43.5 61.75 14.b.2
Perc MaxPer RNDMaxPer AreasList Shape_Leng Shape_Area
1 100.00000000 100.0000000 100 14.b.2 3 0.5
2 84.12674145 84.1267414 84 14.b.2 3 0.5
3 24.99803694 24.9980369 25 14.b.2 3 0.5
4 11.97744244 11.9774424 12 14.b.2 3 0.5
5 3.89717625 3.8971762 4 14.b.2 3 0.5
6 0.28578428 0.2857843 0 14.b.2 3 0.5pts <- SpatialPoints(cbind(pred_grid_tmp3$lon, pred_grid_tmp3$lat),
proj4string = CRS(proj4string(shape)))
pred_grid_tmp4$subdiv <- over(pts, shape)$Area_27
# Rename subdivisions to the more common names and do some more filtering (by sub div and area)
sort(unique(pred_grid_tmp4$subdiv))[1] "3.d.24" "3.d.25" "3.d.26" "3.d.27" "3.d.28.1" "3.d.28.2" "3.d.29"
[8] "3.d.32" pred_grid_tmp4 <- pred_grid_tmp4 |>
mutate(sub_div = factor(subdiv),
sub_div = fct_recode(subdiv,
"24" = "3.d.24",
"25" = "3.d.25",
"26" = "3.d.26",
"27" = "3.d.27",
"28" = "3.d.28.1",
"28" = "3.d.28.2",
"29" = "3.d.29"),
sub_div = as.character(sub_div)) |>
filter(sub_div %in% c("24", "25", "26", "27", "28", 2)) |>
filter(lat > 54 & lat < 59 & lon < 22)
# Add ICES rectangles
pred_grid_tmp4$ices_rect <- mapplots::ices.rect2(lon = pred_grid_tmp4$lon, lat = pred_grid_tmp4$lat)
plot_map +
geom_raster(data = filter(pred_grid_tmp4, year == 1999), aes(X*1000, Y*1000, fill = depth)) +
facet_wrap(~sub_div)Warning: Removed 969 rows containing missing values or values outside the scale range
(`geom_raster()`).
pred_grid_done <- pred_grid_tmp4 |> dplyr::select(-subdiv)Check
#zero valued env variables
zeropgplot <- pred_grid_done |>
pivot_longer(cols=c(depth, oxy, sal, temp), names_to = "variable", values_to = "value" ) |>
mutate(value_sc = scale(value), .by = variable) |>
filter(value < 0)
plot_map +
geom_raster(data = zeropgplot, aes(X*1000, Y*1000, fill = value_sc)) +
#facet_wrap(~year) +
facet_wrap(~variable) +
scale_fill_viridis() # seems plausible except for minus temp out at sea. Warning: Removed 2483 rows containing missing values or values outside the scale range
(`geom_raster()`).
plot_map +
geom_raster(data = zeropgplot |> filter(variable == "temp"), aes(X*1000, Y*1000, fill = value_sc)) +
facet_wrap(~year) +
scale_fill_viridis() # 1996 and 2006 was cold yearsWarning: Raster pixels are placed at uneven horizontal intervals and will be shifted
ℹ Consider using `geom_tile()` instead.Warning: Removed 2155 rows containing missing values or values outside the scale range
(`geom_raster()`).
# cod density
plot_map +
geom_tile(data = pred_grid_done |> filter(between(year, 1993,2021)), aes(X*1000, Y*1000, fill = exp(log_density_cod))) +
scale_fill_viridis(trans = "sqrt")Warning: Removed 28101 rows containing missing values or values outside the scale range
(`geom_tile()`).
Save
saveRDS(pred_grid_done, file = paste0(home, "/data/pred_grid.rds"))